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Abstract 

r-| ' We derive exact dispersion relations for axial and flexural elastic wave mo- 

tion in a rod and a beam under finite deformation. For axial motion we 
consider a simple rod model, and for flexural motion we employ the Euler- 
Bernoulli kinematic hypothesis and consider both a conventional transverse 
motion model and an inextensional planar motion model. The underlying 
formulation uses the Cauchy stress and the Green-Lagrange strain without 
omission of higher order terms. For all models, we consider linear constitu- 
tive relations in order to isolate the effect of finite motion. The proposed 
methodology, however, is applicable to problems that also exhibit material 
nonlinearity. For the rod model, we obtain the exact analytical explicit so- 
lution of the derived finite-deformation dispersion relation, and compare it 
with data obtained via numerical simulation of nonlinear wave propagation 
in a finite rod. For the beam model, we obtain an approximate solution 
by standard root finding. The results allow us to quantify the deviation 
in the dispersion curves when exact large deformation is considered com- 
pared to following the assumption of infinitesimal deformation. We show 
that incorporation of finite deformation following the chosen definitions of 
stress and strain raises the frequency branches for both axial and flexural 
waves and consequently also raises the phase and group velocities above the 
nominal values associated with linear motion. For the beam problem, only 
the inextensional planar motion model provides an accurate description of 
the finite-deformation response for both static deflection and wave disper- 
sion; the conventional transverse motion model fails to do so. Our findings, 
which represent the first derivation of finite-strain, amplitude-dependent dis- 
persion relations for any type of elastic media, draw attention to (1) the 
tangible effect of finite deformation on wave dispersion and consequently on 
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the speed of sound in an elastic medium and (2) the importance of incorpo- 
rating longitudinal-transverse motion coupling in both static and dynamic 
analysis of thin structures subjected to large nonlinear deformation. 



1. Introduction 

The dispersion relation provides a fundamental characterization of the 
nature of wave propagation in an elastic solid. Its derivation for various 
solid configurations, such as beams, plates, surfaces etc., has been key to 
the development of the field of elastic wave propagation, a field that traces 
its beginnings to the seminal memoir of Poisson [1]. In this memoir, and 
in the work of Cauchy [2J, it was revealed that two types of waves exist in 
solids: longitudinal and transverse. Analysis of the dispersive nature of these 
waves in different types of unbounded structures has been a focus of many 
classical studies since then. Brief surveys on the historical development of 
theoretical elastic wave propagation research are provided by Graff [3] and 
Ben-Menahem and a thorough discussion on the reconciliation of theory 
with experiments on mechanical waves is given by Thurston (^j. 

These early investigations as well as the vast majority of contemporary 
studies of wave propagation in elastic solids are based primarily on linear 
analysis, that is, linear constitutive laws and linear strain-displacement rela- 
tionships are assumed (see j^j, j^j, and references therein). The incorporation 
of nonlinear effects has nevertheless been actively pursued and continues to 
attract attention as it allows for a more accurate description of the underly- 
ing motion and facilitates the study of amplitude-dependent wave interaction 
phenomena that do not appear in linear systems (see the monograph on non- 
linear oscillations by Nayfeh and Mook |7| and a recent special journal issue 
edited by Destrade and Saccomandi [8[ for broad listings of references on 
the subject). From an engineering perspective, the effects of nonlinearity on 
the dispersion of waves in waveguides could be utilized to enrich the design 
of materials and structural components. In general, the study of nonlinear 
elastic wave propagation is relevant to nonlinear vibration analysis 



dislocation and crack dynamics analysis [ll|, [12| , geophysical and seismic mo 
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tion analysis 13l.ll4j]. material characterization and nondestructive evaluation 



15j, |16|, biomedical imaging [17J, and others. 

Finite amplitude wave propagation in elastic solids is a subset among the 
broader class of nonlinear wave propagation problems. From a mathematical 
perspective, a formal treatment of finite deformation requires the incorpora- 
tion of an exact nonlinear strain tensor in setting up the governing equations 
of motion 18[. As a result, the emerging analysis should permit large and 
finite strain fields, as opposed to small and infinitesimal strain fields. A large 
portion of research on finite amplitude waves considers initially strained ma- 



terials (see, for example, the early studies by Truesdell [19[ and Green [20 
for extensive discussions on the topic). Among the relatively recent works 
that focused on finite-amplitude plane waves in materials subjected to a 
large static finite deformation include those of Boulanger and Hayes 21], 
Boulanger et al. 22J and Destrade and Saccomandi 23j. Furthermore, anal- 



and Destrade and Saccomandi 
ysis of finite-amplitude waves in solids often involve sm all p arameters or 

for example, 



asymptotic expansions (see Norris [24| for a review). Auld [25 
and subsequent _studies by de Lima and Hamilton [26j], Deng 
vastava et al. 
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treated the problem of finite-strain waves using normal 
mode expansion and forced response calculations. This approach is based on 
perturbation theory and is therefore limited to small amplitudes. Focusing 
on rods, modeled in one dimension or higher, many studies considered finite 
amplitude waves for both incompressible and compressible materials (e.g., 
[5], H^-35|). Of particular relevance is a recent investigation by Zhang and 
Liu [36] in which an exact equation of motion for a rod and an approximate 
equation of motion for an Euler-Bernoulli beam were derived under the con- 
dition of finite deformation. It is rather common in the study of structural 
waveguides to omit (for simplicity) some terms in the strain tensor at the 
outset, and among the strain tensor terms that are retained some emerging 
high order terms associated with the deformation of the cross-sectional frame 
are often neglected. 

A broad overview of the listed references and other finite-deformation- 
based studies in the literature reveals that generally the interest has been 
in obtaining spatial/temporal solutions, or solutions at certain physical lim- 
its, rather than complete dispersion relations that directly relate frequency to 
wavevectors (or wavenumbers) . In this paper, we provide an exact analysis of 
static deflection and elastic wave dispersion in a rod and a beam under finite 
deformation with no modeling limitation on the amplitude of the deflection 
or the travelling wave. Starting with Hamilton's principle we consider both 
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axial deformation (to represent longitudinal motion in a rod) and flexural 
deformation (to represent transverse or longitudinal-transverse motion in a 
beam). Focusing on homogeneous waveguides with constant cross-section 
and non-dissipative, isotropic material properties, we derive the exact equa- 
tion of motion and dispersion relational For the rod problem, we subsequently 
obtain the explicit frequency versus wavenumber solution. For the flexural 
beam problem, we obtain a solution by numerically finding the roots of the 
derived dispersion relation. In all our derivations, we use the Cauchy and 
Green-Lagrange definitions of stress and strain, respectively. Furthermore, 
all terms in the nonlinear strain tensor are retained and no high order terms 
emerging from the differentiation are subsequently neglected. In order to iso- 
late the effect of finite motion, we consider linear constitutive relations. The 
proposed methodolgy, however, is applicable to problems that also exhibit 
material nonliearity. To validate our theoretical approach, we first consider 
static deflection, in both a rod and a beam, and make comparisons with finite 
element (FE) solutions based on exact elasticity theory. We then consider 
the dynamics regime, where for the rod problem we examine wave propaga- 
tion in a finite rod by means of standard finite-strain numerical simulations 
(also using FE analysis) and compare the full-spectrum response with our 
derived dispersion relation. 

The dispersion formulations we develop in this work, which represent 
the first derivation of finite-strain, amplitude-dependent dispersion relations 
for any type of elastic media, allow us to examine (1) the effect of finite 
deformation on the frequency and phase/group velocity dispersion curves in 
rods and beams and (2) the role that longitudinal-transverse motion coupling 
plays in the underlying nonlinear mechanics of flexural beams. 

2. Rod and beam kinematics based on finite deformation 

2.1. Classification of rod and beam kinematics 

In this paper we study geometrically nonlinear wave propagation in a 
rod and in a flexural beam, assuming slender structure for each. The rod 
case is also equivalent to the problem of propagation of plane longitudinal 
waves along a single direction in a medium in which no lateral boundaries 



1 A derivation of finite-strain dispersion for a beam model based on an approximate 



equation of motion was published earlier by the authors [37 1 
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are imposed [i.e., in the linear regime, the model represents a nondispersive 
one-dimensional (ID) medium]. To put the kinematical description of the 
problems we consider in context, we show in Tabled] the categories of a class 
of theories for the treatment of a rod/beam, ranging from 1 to 4 degree-of- 
freedom (DOF) systems. The table lists longitudinal motion for a rod (one 
model), transverse motion for a beam (one model), planar motion for a beam 
(two models) and spatial motion for a beam (two models). Each model is 
described fundamentally by the number and type of variables, the form of 
the exact displacement field, A, and the constraint (if any) that applies to 
the displacement field variables. 

The rod model admits longitudinal motion, tension or compression, with 
u denoting its axial displacements. In the conventional beam model, the 
axial displacements are assumed to be zero. Hence, the only variable is 
v, which denotes the lateral displacements. The planar beam models admit 
both u and v displacements and provide a coupling relationship between these 
two variables. Figure [1] provides graphical descriptions of the kinematics of 
a planar beam showing the variables in the s — y plane, where s is the 
Lagrangian longitudinal coordinate. The third axis perpendicular to s and 
y is denoted by z. In a general planar beam model, the following kinematic 
relationship holds: e = r — 1, where e is the axial strain of the centerline and 
r = ^/(l + u') 2 + v' 2 . The inextensional planar beam model is a special case 
of the general planar beam model, in which an inextensionality constraint is 
applied such that the axial strain of the beam's centerline is assumed to be 
zero. This provides us with the following relationship: 



u 



Vl-v' 2 -l. (1) 



Application of an inextensionality constraint to a beam with any type of 
boundary conditions is known to be generally adequate in the absence of 



large axial forces (Crespo da Silva and Glynn 38]; Nayfeh and Pai |39l 
Lacarbonara and Yabuno [40]). In all planar beam models, it is assumed 
that the cross-section remains plane and perpendicular to the centerline after 
elastic deformation. In addition, it is assumed that in-plane and out-of- 
plane warpings of the cross-section do not occur. These assumptions are 
justified, even for large finite deformation, in light of our focus on beams 
with aspect ratios high enough for the Euler-Bernoulli kinematic hypthesis 
to be valid. In the spatial beam models, the beam admits four degrees of 
freedom, u, v, w, and 7, where w denotes the lateral displacements in a 
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direction perpendicular to u and v, and 7 is the torsion angle of the cross- 
section. Here the displacement field is written in terms of the matrix Rbs, 
which is defined as 



R 



us 



cos a — sin a 
sin a cos a 
1 



cos sin (5 

1 
— sin (3 cos j3 



1 

cos 7 — sin 7 
sin 7 cos 7 



in which a 
defined as 



arctant>'//i and (3 = — arctanw'/r, where h is an agent variable 



h=l+u', (2) 

and r = \/h 2 + v' 2 . We note that a and /3 are dependent variables represent- 
ing the in-plane and out-of-plane deflection angles, where the torsion angle 
7 is an independent variable. In a similar manner to the planar models, an 
inextensionality constraint can be applied to produce the inextensional spa- 
tial beam model. The advantage of the inextensionality constraint is that it 
reduces the number of DOF of a system by one, as indicated in Tabled] The 
cases we consider in this work (rod, conventional beam and inextensional pla- 
nar beam) are all 1 DOF systems. Future work will address finite-strain wave 



dispersion behavior of higher DOF systems (see |41( for a static treatment of 
higher DOF sytems). 



2.2. Finite strain fields of single DOF models 

In Section [3j we derive the equation of motion for our three models on 
the basis of the Green-Lagrange strain tensor. In the following subsections, 
the form of this tensor is given, for a rod and for a beam. Other measures of 
strain may be chosen. However, in this work, we only use the Green-Lagrange 
strain owing to its wide use in the literature. 

2.2.1. Rod: Finite strain field 

Introducing A to generally represent elastic displacement, the Green- 
Lagrange strain field in a rod is given by 
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Table 1: Rod and beam kinematic models 



Description 



Motion/DOF Indep. Variable(s) Displacement field, A 



Constraint 



Rod 

Conventional 
beam 

Inextensional 
planar beam 

Planar beam 

Inextensional 
spatial beam 

Spatial 
beam 



Longitudinal/ 1 
Transverse/ 1 

Planar/ 1 

Planar/ 2 
Spatial/3 

Spatial/4 



■ft 



u, v 



it, v 



u, v, w, 7 



u, v, w, 7 



ft 



,/2 



v-y 



-yv'/VlT 

(1-i/vTT" 



,'2\ 



u 

V 

w 

u 
v 
w 



u — yv' jr 
v — y{r — h)/r 

u — yv' jr 
v — y{r — h)/r 

{Rbs - I) 
(Rbs - 1) 



r - 1 = 



\/r 2 + w' 2 — 1 



No Deformation 



X 



Infinitesimal (Small Amplitude) Deformation 




Finite (Large Amplitude) Deformation 




(a) 



Ay 




(b) 



s 

-t> 



Figure 1: Kinematical representation of a planar flexural beam model: (a) infinitesimal 
versus finite deformation, (b) illustration of the variables u and v and the rotation angle 
a, noting that r = 1 + e where e is the strain of the centerline. 



<9A 1 /dA\ 2 



where the first and second terms on the RHS represent the linear and non- 
linear parts, respectively. As seen in Table [fl the elastic displacement field, 
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A, is scalar for a rod and is equal to the axial displacement, u. 



2.2.2. Beam: Finite strain field 

In analogy to fl3]), the Green-Lagrange strain tensor for a beam is 

_ 1 SdAj | dAj | dA m dA m \ 
13 2 \dxj dxi dxi dxj J ' 

where i and j take the value of 1 or 2, in which x\ = s and x<i = y (see Fig. 
[Hand Tabled]). For conventional and inextensional beams, A represents the 

elastic displacement vector, which is \—yv'/ yl + v 12 v — y(l — 1/yi + v' 2 
and [it — yv'/r v — y(r — h)/r] T , respectively. 



3. Equations of motion based on finite deformation 

In this section we derive the equation of motion for a rod (Section 13. II) . a 
conventional beam (Section 13. 2p and an inextensional planar beam (Section 
13. 3p . In our derivations, we consider non-dissipative, isotropic media and 
assume constant material and geometric properties. Furthermore, we ignore 
the effect of lateral inertia on the longitudinal motion in all models. 

3.1. Rod: Longitudinal motion 

Using Hamilton's principle, we write the equation of motion for a rod 
under uniaxial stress as 

[ {ST - 5U e + 5W nc ) dt = 0, (5) 
Jo 

where T, U e , and W nc denote kinetic energy, elastic potential energy and 
work done by external non-conservative forces and moments, respectively, 
and t denotes time. Using integration by parts, the variation of kinetic 
energy is 



5T = -pA / (uSu)ds. (6) 
Jo 

where / denotes the length of an arbitrary portion of the rod, and p and A 
denote densoty and cross-sectional area, respectively. The variation of elastic 
potential energy is also obtained using integration by parts and is given as 
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5U e = [ [ (a5e)dAds, (7) 

where a and e are the axial stress and axial strain, respectively. In this work 
we base our analysis on the Cauchy stress and consider a linear stress-strain 
relationship as given by Hooke's law, that is, a = Ee, where E is the Young's 
modulus. Using Eq. (CO), and with the aid of integration by parts, we can 
now write the variation of elastic potential energy as 



5U e = J l~EAh{h 2 - Jds, 



(8) 



where v! = du/ds = u s and h is as defined in Eq. (j2J). The variation of the 
work done by non-conservative forces and moments is given in terms of the 
variation of axial deformation, u, and the distributed external axial load, q u , 

SW nc = [ (q u 5u)ds. (9) 



o 



Substitution of Eqs. (JBJ), (jSJ) and OH]) into Eq. (jSJ) produces the equation of 
motion ( TTOj) and the companion boundary conditions (fTT|) : 



jT |jT (A 1 8u)ds+ (Bt5u + BM) 



s = I 
s = 



dt = 0, (10) 



(£i = or u = 0) and (Bi = or u' = 0) . (11) 

We can now write an exact nonlinear equation of motion of a ID rod under 
finite deformation as 



A x = : pAu = ^EA(3h 2 - l)u" + q u . (12) 
The section loads are 



Bi = ]-EAh(h 2 - 1) and B x = 0, (13) 
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where Bi is the axial force. If the axial deformation is infinitesimal, then vl 
is small and from Eq. ([5]), h « 1. Substitution of h = 1 into Eq. (1T2]) leads 
to 

Ai = : pAu = EAu" + q u . (14) 

which is the equation of motion describing infinitesimal axial deformation. 

3.2. Conventional beam: Transverse motion 

Using Hamilton's principle, we state the equation of motion as follows: 

[ {5T -5U 9 -5U e + 5W nc )dt = 0. (15) 
Jo 

We express the elastic deformation using the flexural displacement v as well 
as the angular displacement a. These two elastic coordinates are related 
by the following holonomic constraint: a = tan _1 (dw/ds) = tan -1 i/. The 
variation of kinetic energy, with the aid of integration by parts, is 



5T 



( 1 












s = 




s = 



I 

' 



(16) 



where J is the second moment of cross-sectional area, defined as 



J = j y 2 dA = J z 2 dA = ^a 4 . (17) 

In Eq. f|T6|) . r is defined as r = vl + v ' 2 . The variation of gravitational 
potential energy is 

5U 9 = gpA [ (5v) ds. (18) 
Jo 

where g denotes acceleration due to gravity. As in the rod case, we use the 
Cauchy stress and the standard linear Hooke's law in our formulation. Since 
we are not permitting in-plane and out-of-plane warpings of the cross-section, 
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we set the Poisson's ratio, u, to zero in the constitutive relationship. The 
variation of elastic potential energy can now be written as 

SU e = (cr ss 5e ss + o-yydtyy + 2 a sy 5e sy ) dAds . (19) 

Jo J A 

After substituting Eq. fll]) and the constitutive relationship into ( IT9l) . we 
obtain 

SU e = [ {C vp 8v' + C vpp Sv"}ds. (20) 
Jo 



where the coefficient of Sv ' is 



l„ A/ts « EJ(2-3r 2 )v'(v") 2 EJ f v'(v"Y 
C vp =-EA{v'f + J - f rl { , (21) 



and the coefficient of Sv is 

EJ f {v"f _ EJ(l-3r 2 )v" 

In Eqs. (|2lj) and (|22). J f is the 4 th moment of cross-sectional area, defined 
as 

J f = [ y *dA = [ z A dA = la 6 . (23) 

J A J A 8 

It is also useful to define J c , the higher-order product of cross-sectional area, 



J c = y^VdA = ^a 6 . (24) 
Both Jf and J c have been calculated for a circular cross-section with radius 



a. It can be seen from Eqs. f l23|) and (J23j) that for a circular cross-section 
beam, Jf = 3J C . The variation of the work done by non-conservative forces 
and moments is given in terms of the variation of the flexural displacement 
and the distributed external load, 

SW nc = [ (q v 5v)ds. (25) 
Jo 
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Substitution of Eqs. QSJft, (USD, 02QJ and (USD into Eq. (JTHD generates the 
equation of motion ( l2"B"j) and the companion boundary conditions (T2T1) . 



(A 2 <5t;) ds + (5 2 <5t; + B~ 2 5v') 8 2 q 

V. 



dt = 0, (26) 



(£ 2 = or w = 0) and (B 2 = or v' = 0) . (27) 

Consequently, an exact nonlinear equation of motion for a conventional Euler- 
Bernoulli beam based on finite deformation is 



A 2 = : -pAg - pAv + pj {^j + q v + (C vp )' - (C vpp )" = 0. (28) 

After substituting the C vp and C vpp coefficients form Eqs. (I2TI) and (1221) and 
simplifying, we obtain the closed-form equation of motion for a conventional 
beam, 



A 2 = 0:-p^-pAt; + pJ(a) + g 3 + ^ + 71 7^ 

2 (1 + [v'Y) 

EJ{2 + 3{v') 2 )v^ EJ{1 + A(v') 2 - 9(v') 4 )(v") 3 3EJ f {v") 2 v^ 

2(1 + (i/) 2 ) 2 + 2(l + (i/) 2 ) 4 2(l + (i/) 2 ) 4 

2AEJ f v\v"fv^ 3EJ f (l-9(v') 2 )(v") 5 3EJ f v"(v^) 2 _ 

+ (1 + {v'ff + (1 + (v') 2 f (1 + (vyf ~ ' 

(29) 



The section loads, namely the transverse shear force and the bending mo- 
ment, are respectively 

B * = PJ (^J + Cv P - (C vpp )' and W 2 = C vpp . (30) 
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3.3. Inextensional planar beam: Longitudinal-transverse motion 

Following Hamilton's principle, the inextensionality constraint as stated 
in Eq. (jTJ), and its corresponding Lagrange multiplier, we write 



(ST - 5U g - 5U e + SW nc + 5 (tf.e)) alt = 0. 



(31) 



Accounting for rotary inertia, we write the total kinetic energy of the beam 
as 



T=l/2p f {A(u 2 + v 2 ) + Ja 2 }als. 
Jo 



(32) 



The variation of kinetic energy is then obtained with the aid of integration 
by parts, 



ST = -p[ L {{A[u v]-[-$Ja $Ja]'} 
Jo 

-p[-&J& ] 



Su 
8v 



}als 



Su 
Sv 



s = I 
s = ' 



(33) 



Substituting Eq. fl3j) and the linear constitutive relationship into Eq. (|T9l) 
yields the unknown coefficients of the variation of elastic potential energy in 
Eq. <m. The coefficient of Sv' in Eq. (EDI is 



a 



vp 



J 



Jf(l + 2v' 2 )v" 2 )\ Ev'v" 2 



2(1 -v' 2 ) J (l-v 12 ) 2 



J 



3v 



1/2 



2(1 -v' 2 ) 1-v 



Similarly, the coefficient of Sv" in Eq. (12111) is 



a 



vpp 



Ev' 2 v® 



,12 



(34) 



(35) 



2(1 -v' 2 ) 

Substitution of Eqs. (j35]l . (1TH]) . (120|1 and fl55|> into Eq. (J3TJ generates the 
equation of motion ( 136|) and the companion boundary conditions (137|) . 




(A 2 fo) ds + (B 2 5v + B 2 5v') 



L 



dt = 0, 



(36) 
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(B 2 = or v = 0) and (B 2 = or v' = 0) . (37) 

Consequently, an exact nonlinear equation of motion for an inxetensional 
Euler-Bernoulli beam based on finite deformation is 



A 2 = 



—pAg — pAv + pj 



a 



+ q v + (C vp )' - (C vpp ) /f = 0. (38) 



After substituting the C vp and C vpp coefficients form Eqs. (JMj) and fl35l) and 
simplifying, we obtain the closed-form equation of motion for an inextensional 
planar beam, 



A 2 = : - pAg - pAv + pJ 



a 



,'2 



(pAu - q 2 )ds 



+ Q3 



EJvW 4EJv'v"v^ EJ{l + 3(v') 2 )(v"f 3EJ f (v"] 



2 7,( 4 ) 



■ A2 



A2\ 



\2\ 3 



2(1 

UEJ f v'(v") 3 v® 3EJ f (l + h(v') 2 )(v"f 3EJ f v"{ 



A2\ 



A2\ 



2(1 



A2\ 



A2\ 



= 0. 

(39) 



The transverse shear force and bending moment section loads are found from 
Eqs. (140 ap and (140bj) . respectively. 



Bo 



EJv'iv 



n\2 



1 



A2\ 



EJv® 



(v') 



,12 



(pAu — q 2 )ds 
3EJ f (v") 2 v^ 3EJ f v'(v 



/A 4 



<\2 



2(1 



A2\ 



2(1 



^3' 



(40a) 



B 2 = EJv"+ A ; 



2(1 - w 



(40b) 



It is noticeable from Eqs. ff23|) and f|T7|) . that Jf is quantitatively small 
compared to J when the beam is slender (i.e., when a is small). As a result, 
we can neglect the terms containing Jf in Eqs. ( l39l) . ( 140aP and ( I40bl) (fur- 
ther discussion on this omission is provided in Section 16. 3j) . If these terms 
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are dropped from Eq. (1391) only three elastic terms remain, namely — fzrrk 



l\2 ; 



— iE -, Jv ', v ',lfJ± and — EJ ^ l ^ < flJ}i ^ ■ The resulting equation of motion is rela- 

(i—(v'Y) (l— (i/j ) 

tively simple yet still adequate for large elastic deformation. 



4. Static deflection under finite deformation 

In this section we consider the static deflection problem and hence omit all 
the inertial terms in the rod and beam equations of motion (Eqs. ffl2l) . (T28l) 
and (1311 ). Figure [2]presents the static response of a rod due to the application 
of a force load, q u (Section 13.11 ). We show the response following our finite- 
deformation formulation and that obtained using the COMSOL Multiphysics 
software package [42| which is based on a three-dimensional (3D) nonlinear 
FE model of a rod following the theory of elasticity^] The results are in 
excellent agreement and thus provide a validation to our rod formulation for 
static response. Figure [2] also shows that omission of the Poisson' ratio has 
a negligible effect on the nonlinear static response. The deflection under the 
infinitesimal strain assumption is included for comparison. 

In Fig. [3] we consider a clamped-free beam, with a length of L, under 
a uniform external load, q v . In the figure we show the results obtained by 
our finite-deformation formulations (following the conventional transverse 
motion model and the inextensional planar motion model) and the result 
obtained by the COMSOL Multiphysics software package (based on the same 
3D nonlinear FE model used for the rod). Once again, the deflection under 
the infinitesimal strain assumption is included for comparison. 

It is clear that the inextensional planar beam theory model agrees very 
well with the nonlinear FE solution based on the theory of elasticity. This 
is a significant result considering that the inextensional planar beam theory 
model is only a 1 DOF model. The conventional beam theory model on 
the other hand perfoms poorly and fails to even qualitatively capture the 
basic trend of the response. Here we note that in the conventional beam 
theory model there is a term that includes "EA" (Eq. (T2T|) ). Yet, in the 
inextensional planar beam theory model any term including "EA" drops out 
as a consequence of the derivations. Should we forcefully drop the term that 



2 The COMSOL model used is "solid mechanics with geometric nonlinearity" . A fine 
mesh density (12,000 finite elements) was utilized and the calculations were stationary 
with a maximum of 25 iterations in each step. 
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Displacement, u/L 

Figure 2: Static deflection of a rod. 

features "EA" in Eq. (T2Tj) we find that the response improves significantly 
although still does not match the accuracy of the inextensional planar mo- 
tion model. For the wave propagation problem studied in Sections [5] and El 
omission of the term including "EA" will also improve the performance of 
the conventional transverse motion model although in that case the response 
will not qualitatively follow the solution of the inextensional planar motion 
model (and therefore this omission will not be illustrated in the results we 
present in Section [6]). 

5. Wave dispersion under finite deformation 

In this section we derive exact amplitude-dependent dispersion relations 
for harmonic wave propagation in a rod and in a flexural beam under finite 
deformation. We also obtain the explicit frequency versus wavenumber solu- 
tions - exactly for the rod problem and approximately for the flexural beam 
problem. 
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Position, s/L 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 




(COMSOL nonlinear FE solution) 



Figure 3: Static deflection of a clamped-free flexural beam (q v L 3 /8EJ = —0.5). 

5.1. Rod: Exact analytical dispersion relation and solution 
Using Eq. ([2]), we rewrite Eq. (fT2|) as 

u-clu''=\[^l{u'f + cl{u'f]\ (41) 

where Cq = \J~E~Jp. Equation (HT1) is integrable despite it being nonlinear. 
Differentiating ( I4TT) once with respect to s gives 

(uy-clu^ = l[3cl(ur + cl(ur]''. (42) 

Defining u = u' and r = ut, where u is the frequency of a harmonic wave, 
Eq. (1421 becomes 
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u 2 u TT - c 2 u" = i [3c 2 (u) 2 + c 2 (uf}" . (43) 

Defining z = \k\s + r, where k is the wavenumber of a harmonic wave, we 
rewrite Eq. f|4"3l) as 

u?u zz - c 2 k 2 u zz = ^k 2 [3c 2 (u) 2 + c 2 (u) 3 ] zz , (44) 

where now the explicit dependency on s and r has been eliminated. Inte- 
grating Eq. (pHl) twice leads to 

u 2 u - cln 2 u = 1 -k 2 [3c 2 Q (u) 2 + c 2 {uf] , (45) 

or, 



(co 2 - c 2 y) u - ^ (3m 2 + m 3 ) = 0. (46) 

We note that in our integartion of Eq. f l4"4"|) we get nonzero constants of inte- 
gration in the form of polynomials in z. Since these represent secular terms 
we have set them all equal to zero in light of our interest in the dispersion 
relation. Selecting the positive root of Eq. ( l4"6j) we get 

_ W = zi+yTTgTgg) (47) 

Since u = u s , we recognize that u = \k\u z and therefore Eq. f )4T|) represents 
a first order nonlinear ordinary differential equation with z and u as the 
independent and dependent variables, respectively. 

Now we return to Eq. (I4ip and consider for initial conditions a sinusoidal 
displacement field, with amplitude B and a zero phase in time, and a zero 
velocity field. Following the change of variables we have introduced, these 
initial conditions essentially correspond to the following restrictions on the 
u(z) function given in Eq. ( |47|) : 
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fi(0) = \k\B, u z (0) = 0. (48) 

We note that since z describes a space-time wave phase, the restrictions 
given in Eq. (JJSJ) correspond to an initial wave phase and hence may be 
viewed as "initial conditions" on the wave phase. The importance of these 
initial conditions is that they incorporate the effect of the wave amplitude, 
B. Applying Eq. ( |4"8"j) to Eq. (J4"T|) allows us to use the latter to solve for to 
for a given value of k at z = 0. Thus we obtain an exact dispersion relation 
for wave propagation in a rod under finite deformation, which is 



,2 + W\k\ + (Bk) 2 . . 

u(k; B) = y — i L UiQU (49) 

where U[ n f is the frequency based on infinitesimal deformation, 

Winf(«) = Cq\k\. (50) 

By taking the limit, lim u(k;B), in Eq. (I4"9l we recover Eq. (1501) which is 

B— >o 

the standard linear dispersion relation for a rod. 

5.2. Conventional beam: Exact analytical dispersion relation 

Differentiating Eq. (1281) once with respect to s and omitting the external 
loads and gravity yields 

-pA(v)' + pj (£\ +{C vp )"-(C vpp )"' = 0- (51) 
Defining v = v', r = cot and z = ks + rj§ we re-write Eq. (15T|) as 

- P Au 2 {v zz ) + P Jk 2 u 2 (^f) + K 2 (C vp ) zz - k\C vw ) zzz = 0. (52) 
Integrating Eq, (1321 twice with respect to z gives us 



3 From this point onwards, we replace |k| with k since only even-ordered functions of 
the wavenumber appear in the derived equations of motion. 
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— pAu 2 v + pJk 2 uj 2 ^- + K 2 C vp — K 3 (C V p P ) z = 0, (53) 
where a = tan ?; and 

C vp = l/2EA{vf + EJK 2 {2 ~ 3r ^ z)2 - EJ^-^f, (54a) 

C vpp = EJfJ&g- - EJ K {1 -^ , (54b) 

r = y/l + (v) 2 . (54c) 

In deriving Eq. (|53|) we set the constants of integration to zero and thus 
drop out the secular terms as we did for the rod model. Furthermore, we 
recognize that v = v s and therefore v = \k\v z . In an analogous manner to 
the rod problem (see Sec. 15. li . we write the following initial conditions for 
the wave phase: 

u(0) = Bk, v z (0) = 0. (55) 

Eq. (I53p . which is a function of B, can be solved numerically using a standard 
root finding technique to obtain the explicit dispersion relationship between 
frequency and wavenumber for a conventional beam under finite deformation. 

5. 3. Inextensional planar beam: Exact analytical dispersion relation 

Differentiating Eq. ( )38|) once with respect to s and omitting the external 
loads and gravity yields 

-pAiy)' + Pj(^) +(C vp )"-(C vpp )"' = 0. (56) 
Defining v = v ', r = cut and z = ks + r, Eq. ( 156]) is re-written as 

-pAu 2 (v zz ) + pjK 2 u 2 (^Yl + k2 ( C v P )zz - K 3 (C vpp ) zzz = 0. (57) 
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Integrating Eq. (157j) twice leads to 



-pAw 2 v + pJ/tV^f + k 2 C vp - k 3 (C vpp ) z = 0, (58) 
where a = tan _1 w and 



a 



up 



J 



K 2 J f (l + 2v 2 )v 2 )\ Ek 2 vv 2 



2(1 



V 2 ) 



V 2 ) 2 



J + 



3v 



r,2 



2(1 



V 2 ) 



Ek 2 v 2 v z 
1-v 2 



(59a) 



Ekv z . (59b) 

In deriving Eq. (1581) we set the constants of integration to zero in order 
to drop out the secular terms as we did for the rod and conventional beam 
models. Furthermore, and as we stated for the conventional beam model, 
v = v s and therefore v = \k\v z . Following the same approach as in Sections 
l5.1l and l5T2| we use the initial conditions stated in Eq. fl55|) for the wave phase. 
Subsequently, Eq. fl58|) . which is a function of B, can be solved numerically 
using a standard root finding technique to obtain the explicit dispersion 
relationship between frequency and wavenumber for an inextensional planar 
beam under finite deformation. 



/ JfK 2 V 2 \ 

vpp ~ V 2(1 - v 2 ) ) 



6. Results and discussion 

6.1. Finite- deformation dispersion curves 

Two amplitude-dependent finite-deformation dispersion curves for an in- 
finite rod based on Eq. (jl9|) are plotted in Fig. HI These dispersion curves 
provide an exact fundamental description of how an elastic harmonic wave 
locally, and instantaneously, disperses in an infinite rod under the dynamic 
condition of amplitude-dependent finite deformation. Superimposed in the 
same figure is the dispersion curve based on infinitesimal deformation, i.e., 
Eq. fl50|) . We observe that the deviation between a finite-deformation curve 
and the infinitesimal-deformation curve increases gradually with wavenumer, 
and the effect of the wave amplitude on this deviation appears to be rela- 
tively steady - in the wavenumber range considered - as B is increased. More 
thorough inspection of this deviation is provided in Section 16.21 
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In addition, the results from standard FE simulations of a finite version 
of the rod with length L are presented. A prescribed sinusoidal displacement 
with frequency Cj and amplitude B, i.e. u(L,t) = B sin Cjt, was applied to 
the tip of the rod with free-free boundary conditions. The finite-deformation 
FE model consisted of 60 piecewise linear elements with equal lengths, and 
each node consisted of two degrees of freedom, u and v! . Equal time steps of 
10~ 4 [s] were considered in the numerical integration which was implemented 



using MATLAB's odell3 solver (43j. The wavenumber has been recorded by 
observing the wavelength after one period of temporal oscillation of the tip 
(i.e., excitation point) of the rod, and plotted as a function of frequency Cj for 
two given amplitudes. This recording was possible because at the vicinity of 
the excitation point the wave's harmonic form was effectively still maintained 
during the first oscillation cycle. The data points from this simulation (of 
a finite rod) match very well with the analytically derived exact dispersion 
curve (corresponding to an infinite rod). While the wave considered at the 
tip of the excited rod will evolve, under finite strain, into a complex form 
as it propagates into the rod, this correlation provides a validation that a 
given harmonic wave will locally and instantaneously disperse in a manner 
exactly as described by Eq. (H9|) . A second set of independent simulation 
data is also shown in the figure. This set corresponds to the response due 
to initial sinusoidal displacements at a prescribed wavenumber (applied at a 
state when the rod is at rest). Here we measured the frequency of oscillation 
as the wave propagates considering the first temporal cycle. The simulation 
parameters are the same to those used in the runs where an initial end- 
point displacement is prescribed. As illustrated in the figure, this simulation 
provides yet another validation of the analytical dispersion relation given by 
Eq. (HHD- 

In Fig. we present dispersion curves for the conventional and inex- 
tensional planar beam models. For each of the models, we show results for 
two different amplitudes, B. In the interest of generality, the amplitude is 
considered in relation to the beam's cross-sectional radius, a. The explicit 
relationship between frequency and wavenumber is obtained by numerically 
solving for the roots of Eq. (153 p for the conventional beam model and the 
roots of Eq. f )58|) for the inextensional planar beam model. Also shown in 
Fig. [5] is the dispersion curve based on infinitesimal deformation, which is 
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— Infinitesimal strain - exact dispersion relation 
Finite strain - exact dispersion relation 

- - 5 = 0.1 




Wavenumber, k 

Figure 4: Frequency dispersion curves for a rod. The finite-strain dispersion relation is 
based on Eq. (|49|) : the infinitesimal strain dispersion relation is based on Eq. ([50j) . 



2 

The finite-strain curve of the conventional beam model shows a rapid change 
in slope at low frequencies compared to the infinitismal-strain curve, which 
suggests poor peformance as was seen in the static case. The inextentional 
beam response on the other hand shows asymptotic convergence with the 
linear dispersion as u — > and the deviation grows slowly with k. We can 
also observe the effect of the wave amplitude on the disperison: as B increases 
the dispersion curve rises at an increasing rate. Further comparison between 
the finite-deformation and infinitesimal-deformation responses is provided in 
Section 16.21 
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1 2 3 4 5 

Wavenumber, k (rad/m) 



Figure 5: Frequency dispersion curves for a conventional beam and an inextensional planar 
beam. The finite-strain curves represent exact solutions of Eqs. (|53|) and (|58p. respectively, 
evaluated by numerical root-finding. The infinitesimal-strain curve (based on Eq. (|60p ) is 
also shown. 

Figure |6] shows the group velocity versus wavenumber curves for the same 
cases considered in Figs. E2 and El The group velocity is defined as c g = 
du/dn. The figure puts on view the significant qualitative transformation 
that emerges when finite deformation is incorporated in wave propagation 
analysis. For axial waves the group velocity, which is otherwise constant, 
is now varying with k. And for flexural waves the group velocity follows a 
nonlinear function with k in contrast to the linear trend for the infinitesimal- 
deformation model. Needless to say, the intensity of these finite-deformation 
effects on the group velocity increases with B as mentioned in our discussion 
of the frequency dispersion curves. A phase velocity dispersion diagram may 
also be generated, and from Figs. H] and [5] it is clear that the trends will be 
similar to Fig. O 

The dispersion characteristics presented provide an elucidation of the non- 
linear, finite-deformation dynamics of rods and beams particularly for long 
structures or when the focus is on high frequency /short wavelength behav- 



24 



10 





6 -! 












4 -j 


o 
_o 




53 
> 


2 [ 


o 


* 


a 




-2 - 



Conventional 



I beam theory 



Inextensional planar 
beam theory ^ ' 




Infinitesimal strain 

Finite strain (rod, B = 0.1; beam, B = 10a) 
Finite strain (rod, B = 0.2; beam, B = 20a) 



0.5 1 1.5 2 

Wavenumber, k (rad/m) 



2.5 



Figure 6: Group velocity dispersion curves for a rod, a conventional beam and an inex- 
tensional planar beam under finite strain. 



ior where a standard analysis based on superposition of standing waves is 
generally not effective 44]. From a conceptual point of view, an increase in 
group velocities as seen in Fig. [6] is particularly significant in the study of 
waves where the rate and intensity of energy propagation is consequential, 
such as in seismic waves [13|, [14J. It is also relevant to the study of the speed 
of propagation of dislocations and cracks considering that large amplitudes 
have been observed near slip planes and crack tips [111 . |l2j (recalling that the 
phase and group velocities change with deformation amplitude). Future work 
will introduce appropriate extensions in the implementation of our method- 
ology to address the specific geometric and material characteristics of these 
applications. 

From an engineering design perspective, the amplitude- dependent disper- 
sion behavior we observe in Figs. HI |5] and |6] potentially could be of benefit 
to numerous applications; for example, the frequency shifts observed suggest 
possible utilization of regular rods and beams as amplitude-dependent wave 
propagation filters. Further utilization is possible when these properties are 



25 



considered in conjunction with other avenues for design such as the intro- 



duction of periodicity |45H47j . Here we note that the derived finite-strain 
dispersion relations can be incorporated directly into transfer matrix models 
of ID periodic media, which will be the focus of future research. 

6.2. Comparison with infinitesimal-deformation dispersion curves 

Figure [7] presents the percentage deviation in the dispersion curves when 
finite deformation is incorporated compared to when infinitesimal deforma- 
tion is assumed (evaluated at k — ir). Clearly the higher the wave amplitude, 
the higher the deviation. For a rod, the deviation follows a concave curve 
(increases at a decreasing rate); whereas for an inextensional planar beam it 
follows a convex curve (increases at an increasing rate). The level of deviation 
for a given wave amplitude is considerably higher for longitudinal waves in 
a rod compared to longitudinal-transverse waves in the inextensional planar 
beam. 



Rod wave amplitude, B 

0.2 0.4 0.6 0.8 1.0 1.2 1.4 




2.0 



2 4 6 8 10 12 14 16 18 20 

Beam wave amplitude to cross-sectional radius ratio, B/a 



Figure 7: Percentage deviation in frequency when the dispersion curves are based on 
infinitesimal strain versus finite strain (evaluated at k = it). 
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6. 3. Effect of the 4 moment of area, Jf 

In this section we examine the effect of the 4 th moment of cross-sectional 
area, Jf, on the finite-deformation dispersion relations for flexural beams. 
If, by way of example, we consider a beam with a circular cross section 
of radius a, then the terms tq and r/ would each be a function of a, i.e., 
r = a/ J/ A = a/2 and r/ = a/ Jj/A = a 2 /2\/2. Hence we get r//r = a/y/2 
from which we can deduce that any term that contains Jf diminishes when 
the cross-sectional area is small. Fig. [8] shows the effect of omission of 
all terms that contain Jf in both the conventional and inextensional planar 
beam models. From the figure we note that for a ratio of cross-sectional 
radius to wave amplitude, a/B, smaller than approximately 0.025, the effect 
of omitting the terms that contain Jf in Eqs. ( )29l) and ( 1391) is negligible. 
However, for higher values of this ratio, the error introduced by this omission 
increases dramatically. 

7. Conclusions 

We have derived and solved an exact dispersion relation for finite-strain 
wave propagation in a rod and an Euler-Bernoulli beam. This represents 
the first derivation of amplitude-dependent dispersion relations for elastic 
media under finite deformation, providing an analogy to the derivation of 
anharmonic (nonlinear) dispersion relations of crystals in condensed matter 
physics. 

For axial motion (i.e., the rod problem), the closed- form exact frequency 
versus wavenumber solution was derived as a function of wave amplitude. 
The solution matches very well with data points obtained from high-resolution 
FE simulations of a finite rod under finite deformation. For flexural motion 
(i.e., the conventional and inextensional planar beam problems), an approxi- 
mate solution was obtained by numerically solving for the roots of the derived 
exact dispersion relation. A key aspect of the derivations that lead to the 
dispersion relation for each of these problems involves the introduction of 
a space-time wave phase, z, and the development of a strain relationship 
which at z — 0, with the appropriate initial conditions, gives the amplitude- 
dependent dispersion relation. 

The results show that finite deformation causes the frequency (and there- 
fore the phase and group velocity) dispersion curves to shift for both axial 
and flexural waves in comparison to when infinitesimal deformation is as- 
sumed. The level at which these nonlinear effects take place is dependent on 
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Figure 8: Percentage error in frequency due to omission of terms including Jf in Eqs. (|28[) 
and (|39l) for conventional and inextensional planar beams (evaluated at k — tt). 



the amplitude of the travelling wave, and the degree of shift in the disperi- 
son curves increases with wavenumber. Qualitatively, finite deformation is 
shown to introduce dispersion to an otherwise nondispersive medium. Both 
the rod and inextensional beam solutions display asymptotic convergence to 
the corresponding linear, infinitesimal-strain solution at low frequencies and 
at low wave amplitudes. For our choice of Cauchy stress and Green-Lagrange 
strain, all observed shifts are positive, i.e., the curves rise with finite defor- 
mation. Future work will examine the response based on other stress and 
strain measures, which in itself will serve as a rigorous characterization of 
the spectral properties of these measures. The results also show, for wave 
dispersion as well as static deflection, that a coupling between the longitu- 
dinal and transverse motions is necessary for an accurate description of the 
response. 
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Our analysis provide a means for elucidating the nonlinear dynamics of 
rods and flexural beams particularly for long structures or when the focus 
is on high frequency/ short wavelength behavior. It also has implications on 
the study of wave phenomena in a broader range of problems where large 
elastic motion takes place. Examples where large deformations have been 
experimentally observed include waves produced by cracks and dislocations 
(problems of great importance in seismology, fracture mechanics and plas- 
ticity). In these problems, the standard dispersion relations based on in- 
finitesimal motion cannot predict any alterations in the speed of sound due 
to the intensity of the motion, and hence the need for a finite-strain treat- 
ment as we observed in Fig. O for the ID models considered. While our 
present formulations consider linear constitutive relations (in order to isolate 
the effect of finite motion), the proposed methodology is applicable to prob- 
lems that also exhibit material nonlinearity - which would be necesary for 
the applications mentioned above. From an engineering design perspective, 
the amplitude-dependent shifts in frequency, phase and group velocity in the 
dispersion spectrum provide promising avenues for the development of novel 
materials, structures and devices. 

The presented analytical formalism for dispersion analysis under finite 
deformation lays the foundation for systematic extensions to more complex 
ID models (e.g., rods/beams incorporating lateral inertia in longitudinal mo- 
tion, beams incorporating corotational motion, higher DOF beams) as well as 
multi-dimensional models (e.g., thin plates and shells, 2D plain stress/strain 
models, semi- infinite surfaces, 3D rods and plates, bulk media). 
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